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A discrete model for the sintering of polydisperse, inhomogeneous arrays of cylinders is presented 
with empirical contact force-laws, taking into account plastic deformations, cohesion, temperature 
dependence (melting), and long-time effects. Samples are prepared under constant isotropic load, 
and are sintered for different sintering times. Increasing both external load and sintering time leads 
to a stronger, stiffer sample after cooling down. The material behavior is interpreted from both 
microscopic and macroscopic points of view. 

Among the interesting results is the observation, that the coordination number, even though it has 
the tendency to increase, sometimes slightly decreases, whereas the density continuously increases 
during sintering - this is interpreted as an indicator of reorganization effects in the packing. Another 
result of this study is the finding, that strongly attractive contacts occur during cool-down of the 
sample and leave a sintered block of material with almost equally strong attractive and repulsive 
contact forces. 
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no assumptions about either geometry, topology, homo- 
geneoity, or isotropy of the powder packing are involved. 
However, we have to assume certain contact force-laws 
and, furthermore, restrict ourselves to two dimensions. 

In order to test the approach, a model system in a 
two-dimensional box filled with cohesive-frictional disks 
of different sizes, see Sec. ||, is examined by means of a 
"microscopic" discrete element method (DEM). The mi- 
croscopic interaction model for plastic deformations, fric- 
tion, and cohesion is discussed and time-, temperature-, 
and history-dependent behavior is introduced and ratio- 
nalized in Sec. lU. The long time sintering of a block 



of material is simulated and the results are discussed in 



Sec. [V 



time 



FIG. 1: (Left) Preparation of the sample with pressure on 
both sides. (Right) After the system is relaxed the evolution 
in time of the sample with time is examined. 



II. MODEL SYSTEM 

One possibility to obtain information about the mate- 
rial behavior is to perform elementary tests in the labo- 
ratory. However, because it is difficult to observe what 
is going on inside the material, an alternative way is 
the simulations with the discrete element model (DEM) 
I, The numerical "experiment" 

chosen is a bi-axial box set-up, see Fig. |l|, where the 
left and bottom walls are fixed, and a stress- or strain- 
controlled deformation is applied to the other walls. In 
the future, quantitative agreement between experiment 
and simulation has to be achieved, however, this issue is 
far from the scope of this paper. 

Stress control means that the wall is subject to a prede- 
fined external pressure that, in equilibrium, is cancelled 
by the stress, which the material exerts on the wall. In 
a typical "experiment" , the particles are rapidly com- 
pressed with constant, isotropic pressure and then left 
alone for some time and at constant temperature, typi- 
cally close below the melting point. 

The stress-controlled motion of the walls is described 
by the differential equation 

m^x{t) = F^{t) - p^z{t) - 7wi(i) , (1) 

where is the mass of the wall. Values of large as 
compared to the sample mass lead to slow adaption and 
vibrations, whereas small values allow for a rapid adap- 
tion to the actual situation. Three forces are active: (i) 
the force F^{t) due to the bulk material, (ii) the force 
—p^z{t) due to the external pressure, and (iii) a viscous 
force which damps the motion of the wall so that oscilla- 
tions are reduced. 

Inside the system, N disks with radii a; (i = 1, . . ., N) 
and height h are placed. The radii are drawn from a ho- 
mogeneous distribution with mean oq and relative width 
wq so that Oi/ao S [1 — 1 + w^]. The particle-particle 
interactions and the parameters involved are discussed in 
the next section. 



III. DISCRETE PARTICLE MODEL 

The elementary units of granular materials are meso- 
scopic grains which deform under stress, possibly yield, 
change their properties with time, and can behave differ- 
ent for different temperatures. Since the realistic mod- 
eling of the deformations of the particles is much too 
complicated to allow for a subsequent many-particle sim- 
ulation, we relate the interaction force to the overlap 5 
of two particles, see Fig. |[ In the absence of long-range 
forces, an interaction takes place only if particles are in 
contact and thus (5 > 0. In that case, the forces are split 
into a normal and a tangential component denoted by n 
and t, respectively. 

If all forces /j acting on the particle i, either from 
other particles, from boundaries or from external forces, 
are known, the problem is reduced to the integration of 
Newton's equations of motion for the translational and 
rotational degrees of freedom 

d2 d2 

with the mass of particle i, its position r.; the total 
force /j = ^cfi ' moment of inertia /i, its angular 
velocity uji = d(^j/dt, the total torque ti = J^c^i ^ fi' 
and the center-contact vector l^. The integration of the 
equations of motion is performed with a standard molecu- 
lar dynamics Verlet algorithm together with Verlet-table 
neighborhood search . 

A. Normal Contact Model 

Two particles i and j at positions and rj , with radii 
ai and aj , interact only if they are in contact so that their 
overlap 

S = (a, + ttj) - {vi - rj) ■ n (3) 

is positive, S > 0, with the unit vector n = n.y = (r^ — 
rj)/\ri — rj \ pointing from j to i. The force on particle i, 
from particle j can be written as /j.- = fl^^n+ffA, with n 
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perpendicular to t. In this subsection, the normal forces 
are discussed. 



1. Short time contact model 

As first step, we discuss the time- and temperature- 
independent behavior of the contact forces between a pair 
of particles. For this, we modify and extend the linear 
hysteretic spring model |2^, |3| . It is the simplest ver- 
sion of some more complicated nonlinear-hysteretic force 
laws ^ , which reflects the fact that at the con- 

tact point, plastic (permanent) deformations may take 
place. The repulsive (hysteretic) force can be written as 



fij — 



with ki < k2, see Fig. ||. 




So) 



loading, 

un/reloading, 

unloading, 



(4) 




FIG. 2: (Left) Two particle contact with overlap S. (Right) 
Force law for two springs with stiffness ki and k2 for initial 
loading and subsequent un/reloading, respectively. Attractive 
forces are possible due to the cohesion strength kc- 

During the initial loading the force increases linearly 
with the overlap S, until the maximum overlap iSmax is 
reached, which is kept in memory as a history parameter. 
The line with slope fci thus defines the maximum force 
possible for a given S. During unloading the force drops, 
on the line with slope ^2, from its value at (5max down to 
zero at the force-free overlap 

,5q = (1 - fci/fc2)^inax ■ 

Reloading at any instant leads to an increase of the force 
along this line, until the maximum force is reached; for 
still increasing 6, the force follows again the line with 
slope fci and 5max has to be adjusted accordingly. Un- 
loading below So leads to negative, attractive forces until 
at the overlap 

ka-ki ^ 



the minimum force /min = —kcSmim i-S- the maximum 
attractive force, is obtained as a function of the model 
parameters ki, k2, kc, and the history parameter Jmax. 
Further unloading leads to attractive forces on the branch 

/ = -M. © 

The cone formed by the lines with slope fci and —kc 
defines the range of possible force values. If a force would 
fall outside the cone, it is forced to remain on the limit 
lines. Departure from these lines into the cone takes place 
in the case of unloading and reloading, respectively. Be- 
tween these two extremes, unloading and reloading follow 
the same line with slope k2. Possible equilibrium states 
are indicated as circles in Fig. |2| where the upper and 
lower circle correspond to a pre-stressed and stress-free 
state, respectively. 



2. Viscous dissipation 

In the case of collisions of particles and for large de- 
formations, dissipation takes place due to the hysteretic 
nature of the force-law. For small displacements around 
some equilibrium state, the model does not contain dis- 
sipation. Therefore, in order to allow for stronger dis- 
sipation and thus faster relaxation, a viscous, velocity 
dependent dissipative force in normal direction. 



fi, = loS , 



(5) 



is assumed with some damping coefficient 70. The half- 
period of a vibration around the equilibrium position, see 
Fig. H, can be computed for arbitrary values of fci and 
as long as the overlap fulfills the condition (Jmin < 5 < 
f^max- In this case, fc2 determines the stiffness and one 
obtains a typical response time on the contact level [E6[ , 



tc — — , with Lo = 

UJ 



k2 
mi2 



(6) 



-'mm — , , , ^max ; 
k2 + kc 



the eigenfrequency of the contact, the rescaled damping 
coefficient 770 — 7o/(2"ii2), and the reduced mass mi2 = 
r7iiTO2/(TOi -l- 7712). The time-step of the simulation <md 
has to be chosen such that 

*MD ~ <c/50 

for a proper integration of the equations of motion, so 
that we chose a typical time-step tMD = 7r/(50a;) after 
the model parameters k2 and 70 are specified pl[|. 



3. Stiffness increase with contact area 

In order to account for the fact that a larger contact 
surface leads to a larger contact stiffness, the coefficient 
k2 is made dependent on the maximum overlap history 
parameter <5„iax (and thus on the force-free overlap Sq), 
as long as the overlap is below the threshold 5^™'^ that 
corresponds to the "complete melting" of the particles. 
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Complete melting is here the limit of an incompressible 
liquid that is contained in the model, however, neither 
discussed in detail nor verified for reasons of brevity. 

The overlap S^™^ corresponding to the stress-free fluid, 
is computed such that the volume fraction in the system 
equals unity. The volume fraction of a dense packing of 
rigid spheres is 

J^soiid = Nirhal/VsoUd = 7r/(2\/3) . 
If all the material would be melted, the volume fraction 



VRuid = Nnhal/VRuid = 1 , 
which leads to the ratio of fluid and solid volumes 

Vfluid/Kioiid = flfluid/ao = 7r/(2V3) . 

The minimum radius for the incompressible melt is 
afluid = ao - f^"""* = ao\/Vfluid/"l4oiid ~ 0.9523 ao which 
corresponds to the maximum overlap 



! / 
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fluid 



0.0477 ao 



in two dimensions. In the following, we use (5^""^ = 0.20 
in order to magnify the range of possible non-fluid over- 
laps. 

The stiffness is maximal in the fluid limit for Sq — S^™'^, 
which corresponds to S^ax — Sf^iix ^ fc2'^**""^/(fe — fci), 
and varies between ki and k2 for smaller overlaps, so that 



k2 

ki + (fc2 - fci) 



if 6, 
if S, 



max ^ ^] 



fluid 
max 
fluid 
max 



(7) 

For large overlaps (in the fluid regime), the stiffness and 
the force is thus only dependent on fc2, independent of 
hi. For smaller overlaps both ki and k2 affect the force 
together with the history of this contact. 

The stiffness in the incompressible fluid should diverge, 
however, for reasons of numerical stability, we have to 
limit the maximum stiffness to k2 . Larger stiffness values 
would require smaller time-steps which would reduce the 
simulation efficiency. Therefore, the model does not re- 
ally take the incompressibility of the liquid into account, 
see Fig. |[ In the fluid no plastic deformation can take 
place so that i5o is flxed and cannot be shifted to larger 
force free overlaps ~ in contrast to the hysteretic model 
described above. 

In summary, the hysteretic stiffness model takes into 
account an increasing stiffness with increasing overlap. 
The first loading is plastic with low stiffness, and sub- 
sequent un- and reloading are stiffer because the mate- 
rial was initially compressed. As a consequence, also the 
maximum cohesive force depends on the maximum com- 
pression which was experienced by the contact during its 
history. If the material is compressed so strong that the 
liquid density is reached, the force-free overlap is equal 
to the fluid equivalent overlap and the material behaves 
like a fluid. 




FIG. 3: Force law with varying stiffness fc2(iSmax), according to 
Eq. (^. If So becomes larger than 5*^""^, the stiffness remains 
equal to k2 and the force remains on the corresponding liquid 
branch with slope k2 (dotted line). The cohesion strength 
is maximal for the max imum contact strength and decreases 
with k2, see subsection [HE for details. 



B. Density Temperature Dependence 

If a solid or a liquid (we assume a simple material here 
- not like water - where the density dependence on the 
temperature is continuous through the phase transition) 
is heated, in general, its volume increases so that its den- 
sity decreases. Therefore, we assume a temperature de- 
pendent density of the single particles (disks with radius 
a and height h): 



P(T) 



na'^h 



P(rmolt) + SpT (Tinelt " T) , (8) 



with the density change per unit temperature Spx- This 
corresponds (in linear approximation) to a change of the 
particle radius 



a{T) = a(T,neit)[l - Sut (T^dt - T)] 



(9) 



with the relative change of the radius per unit tempera- 
ture Sax- This approximation can be used if the range of 
temperatures is rather narrow and the changes per unit 
temperature are very small. 

In the following, we use Sut — lO^'^K^^, so that the 
particle radius is changed by 0.01 per-cent if the temper- 
ature is changed by one Kelvin. In the interesting range 
of temperatures between a low temperature (80'^C) and 
some melting point (120°C), the radius changes by one 
per-cent, accordingly. 



C. Contact Temperature Dependence 

For the temperature dependence, we focus on an in- 
homogeneous material with a melting temperature Tmoit 
and assume that the material behaves static, as described 
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above, if the temperature T is much smaller than the 
melting temperature. The behavior of the stiffness ki is 
schematically shown in Fig. ^ as a function of the tem- 
perature. 




T 



FIG. 4; Schematic plot of the stiffness fci as a function of the 
temperature. 



change of fci(T), see Eq. (|T|). In a pre-stressed situation, 
corresponding to a finite confining force at the contact, 
also is shifted in order to balance the confining force 
- but only after fcifJ^ax became smaller than the confining 
force. 




The material becomes softer when T,„eit — T becomes 
small and will lose all stiffness in the limit Tmcit — T <^ 
0. The temperature range in which the melting takes 
place is quantified by Tvar- In the transition regime 
|2moit — r| « TVar, the particles are significantly softer 
than in the cold limit T^dt — T ^ T^ar- In the hot and 
the liquid regime, T — Tmcit ^ Tvar, one has fci ^ and 
the particles lose their nature, however, the 'incompress- 
ibility' is accounted for with a stiffness ^2, as defined in 
Eq. (0), and (5o = 5""''* is fixed. 



1. Increasing temperature 

When the temperature is increased to a rather large 
value, close to the melting point, two particles under 
stress and in equilibrium due to compressive forces will 
lose stiffness and thus will deform stronger so that their 
overlap becomes larger. Therefore, we assume for the 
stiffness coefficient 



fci(T) 



1 + tanh 



y 

^ [1 + tanh(T)] 



T 



(10) 



with the dimensionless temperature difference r, and the 
typical range of considerable temperature dependency 
Tvar, where the temperature has a strong effect on the 
stiffness of the particles, see Fig. Q If experimentally 
available in the future, the function tanh(x) can be re- 
placed by any other function f{x) that decays from unity 
to zero at a; « 0. 

When fci is reduced due to an increase in temper- 
ature we assume that J^^x remains constant, so 
that one obtains a larger force- free overlap 5q{T) = 
[1 — fci (2^)7^2 (<5max)]^max- Thus the material volume 
shrinks due to sintering at the contact level. 

Note that ^2 is not changed directly when kx{T) is 
decreased, see the left panel in Fig. | or Eq. (0). The 
cohesion in this model, however, is directly affected by a 



FIG. 5: Force laws for varying stiffness k\, according to Eq. 
(p^. (Left) If the temperature is increased, fci is reduced 
while 5max remains constant (dashed line, stress-free case). 
(Right) If the temperature is subsequently decreased, fci is 
increased while Sq remains constant (solid line with slope 
fo). 

If, as an example, the material has a melting point 
^meit = 120"C with a range of softening of Tvar = 10°C, 
for a temperature of T = 118°C, the stiffness fci is re- 
duced to 0.6 of its cold-limit value, for a temperature of 
122°C the stiffness is 0.4, and for a temperature of 130°C 
the stiffness is only 0.12. For a temperature of T = IGO^C 
only a stiffness of 3.10"^ remains. The fact that there is 
a remaining stiffness above the melting temperature can 
be attributed to the inhomogeneity of the material, i.e. 
not all the material melts at the same temperature. 



2. Decreasing temperature 

If later in time the temperature is decreased again, 
ki{T) is adjusted according to Eq. (p^, but since the 
melted (sintered) area around the contact point will 
not return to its previous state, we now assume that 
5^ ~ const., so that the maximum overlap increases 
to the value i5~ax(^) = ~ ki{T)/k2\j see the 

right panel in Fig. ^. Therefore, a temperature cycle 
involving a temperature close to the melting temper- 
ature leads to a contact situation similar to the one 
obtained through some larger maximum compressive 
force. The contact deformation/area and thus the force- 
free overlap become larger due to the partial melting of 
the surface and also the stiffness is increased accordingly. 

The temperature dependence thus can lead to changes 
of the stiffness and increases the overlap (deformation) 
of the particles. It is effective only if the temperature is 
changed - the time dependence will be introduced in the 
next subsection. 
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D. Temperature dependence with time 

In the material there are several time dependent pro- 
cesses taking place. Since we are interested in the long 
time behavior of the material, we assume that heat con- 
duction and equilibration take place instantaneously, as 
long as temperature changes are small and slow. The 
realistic simulation of heat-conduction in the sample is 
far from the scope of this study. Therefore, the particle 
size is adjusted directly with the temperature according 
to Eq. d). 

In addition to the direct effect of a temperature change 
on the particle size and the stiffness, the material may 
change its internal, atomistic structure such that defects 
heal and disappear. This effect will occur mainly in the 
regime of high temperatures close to the melting point. 
Therefore, in order to account for such slow microscopic 
processes, a time dependence is introduced that leads to a 
change of the material stiffness fci with time. The change 
takes place extremely slowly with an algebraic time de- 
pendence, so that fci(T, t) lags behind when varying from 
its actual value to the desired, final value ki{T), as de- 
fined in the above equation ([lO|). 

When the temperature is increased from a small value 
To to T, then fci(To) changes to ki(T) following the law 



ki{T,t) ^ ki{T) 



1 



1 



l-fci(To)/fci(T) to 

which corresponds to the rate of change 

dki{T,t) _ \k,{T) ~ k,{T,t)]^ 



(11) 



dt 



= ± 



fci(r) to 



(12) 



with the time scale to on which a typical change takes 
place. Note that ki, ki{T), and ki{T,t) are different, in 
general, and correspond to the maximum fci, the tem- 
perature dependent fci (T) < fci , and the time dependent 
fci(r,t) that tends towards fci(T). The sign in Eq. (|l|) 
is chosen according to the sign of [fci(T) — fci(r, t)] in Eq. 

(0) ii- 

Assume Tq = 20^0 and T = 118°C, so that fci (To) = 1 
corresponds to fci(r) = 0.6. The stiffness as a function 
of time is plotted in Fig. ^for the time constant to = 10s. 

Note that, due to the factor fci(T) in the denominator 
of Eq. (|l2|), the change of stiffness is faster for higher 
temperatures. In the hot limit, changes take place very 
rapidly, whereas in the cold limit changes are extremely 
slow. 

Finally, we note that the adaption/relaxation of 
fci (T, t) to the desired fci (T) value follows the above equa- 
tions in all cases except when the temperature is going 
down and the fci(T) > ki{T,t). In that situation the 
contacts freeze rapidly and thus have to become strong 
as fast as the system cools down. 
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FIG. 6: Variation of the stiffness fci with time, see Eq. (^l|). 
The inset shows the same function with a logarithmic time 
axis - only after about 10'^ s the final stiffness is reached. 



E. Cohesion dependence on stiffness and friction 

The cohesive properties of a particle contact depend on 
the temperature, in so far that a melted contact should 
have weak tensile and compressive strength. Therefore, 
we couple the cohesive parameter kc to the magnitude of 
fci (T, t) , see Eq. (^Tj) , which decreases with temperature 
increasing. In addition, in order to take into account a 
reduced tensile strength of a soft contact with weak de- 
formation and thus small overlap, the cohesion is directly 
related to the stiffness fc2(^max), see Eq. (0). Thus, we 
propose 



kc{Tj t, ^max) 



fcl(T,t) fc2((5niax) 
fci fc2 



kf: 



(13) 



This is an arbitrary choice for the cohesive force factor 
fcc, but as long as no detailed experimental results are 
available, we stick to this empirical law, see Fig. ^ for a 
schematic picture. 



F. Tangential Contact Model 

The force in tangential direction is implemented in the 
spirit of Cundall and Strack |l^] who introduced a tan- 
gential spring in order to account for static friction. Var- 
ious authors have used this idea and numerous variants 
were implemented, see for a summary and discussion. 

Since we combine cohesion and friction and intro- 
duce time and temperature dependencies, it is conve- 
nient to repeat the model and define the implementation. 
The tangential force is coupled to the normal force via 
Coulomb's law, i.e. /* < /i/", where for the limit case 
one has sliding friction and for the case of small forces 
/*, one has static Coulomb friction. The latter situation 
requires an elastic spring in order to allow for a restoring 
force and a non-zero remaining tangential force in static 
equilibrium due to activated Coulomb friction. 
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As a consequence of tlie coliesion force in normal di- 
rection, attractive forces are possible so that Coulomb's 
law has to be modified 

/* < Kr - /min) , (14) 

with the minimum (maximum attractive) force 

r fc2(^max) - fcl(T,0 ^ 

/mm— 1 I 7 /r \ 11 /rn 4- s: \ "max ■ l-LOJ 
1 + fc2(()max)/Kc(j Omax ) 

The equal '—' in Eq. ( p^ ) corresponds to (fully) acti- 
vated, sliding friction, while the smaller '<' corresponds 
to static friction. Thus the tangential force is related to 
the normal force relative to the point of cohesion-failure. 
For normal forces larger than /min, friction is always ac- 
tive and the amplitude is proportional to /" — /min and 
/i. This can lead to a stable equilibrium of the solid also 
at /" « and activated tangential static friction. Note 
that /min tends towards zero for vanishing overlap, so 
that the difference from the original model also vanishes 
for vanishing overlap. 

If /" — /min > 0, the tangential force is active, and we 
project the tangential spring ^ into the actual tangential 
plane by subtracting a possible (small) normal compo- 
nent. @. 

S. = i-n{n-i). (16) 

This action is relevant for an already existing spring, if 
the spring is new, the tangential force is zero anyway, 
however, the change of the spring-length defined below, 
is well defined. Next, the tangential velocity is computed 

vt = Vij - n{n ■ v.ij) , (17) 

with the total relative velocity 

Vij = Vi — Vj + aiTi X LJi + ajti x u3j . (18) 

In the next step we calculate the tangential test-force as 
the sum of the tangential spring and a tangential viscous 
force (similar to the normal viscous force) 

ft^~kti~ltvt . (19) 

with the tangential spring stiffness kt = afc2('^max), with 
a typical stiffness ratio a = 0.2, see and the tan- 
gential dissipation parameter 7^. If |/°| < fc, with 
fc = f^if"' — /min), one has static friction and, on the 
other hand, if |/°| > fc, sliding friction is active. In the 
former case, the tangential test force is incremented 

I' = I + VtAtMD , (20) 

to be used in the next iteration in Eq. (p^), and the 
force /* ~ f° from Eq. ( p^ is used. In the latter case, 
the tangential spring is adjusted to a length which is 
consistent with Coulomb's condition 

4' = -^(/c*), (21) 
kt 



with the tangential unit vector, t = f°/\f°\, defined by 
the tangential spring, and the Coulomb force is used. 
Inserting into Eq. ( p9| ) leads to /* w f^t in the next 
iteration |34| . In short notation this reads 

f ^+mm{fc,\ft\)t ■ (22) 

Note that the tangential force described above is iden- 
tical to the classical Cundall-Strack spring only in the 
limit 7t = and kc — 0. Besides the combination of the 
cohesive and the frictional force, also the tangential dis- 
sipation is non-standard. Furthermore, we remark that 
the cohesion could also be coupled to friction in the sense 
that a broken contact looses its tensile strength when it 
is assumed brittle, so that kc = 0, (if sliding), i.e. if one 
has a sliding contact with /* = m(/" ~ /min) ]35[ |. 

G. Temperature dependence in tangential direction 

In parallel to the change of normal stiffness, the tan- 
gential stiffness is always kept in a constant ratio to k2 
so that 

kt = Q;/C2((5max) , (23) 

since the stiffness in tangential direction is based on the 
same arguments as the material stiffness in normal direc- 
tion. 

The friction is coupled to the temperature dependent 
value of the stiffness fci (T, t) , because friction should not 
be present in a liquid at large enough temperatures, so 
that 

Mr,i) = ^i^M. (24) 
ki 

Thus friction is modified together with the changes in 
normal direction. No new ideas are introduced for the 
tangential forces. 

IV. RESULTS 

In this section, the sintering model is applied to the 
sintering process of a particulate material sample. The 
material is initially a loose powder and first has to be 
prepared at low temperature from time to to time thcat, 
see Fig. |^. The preparation takes place with a system 
as described in section ^ with isotropic external pressure 
P ■= Pw = Px = Pz = 10 or 100. Pressure is measured in 
units of Nm^^ due to the two-dimensional nature of the 
model, i.e. p has the same units as the spring stiffness 
k2 ■ The correct units of the pressure can be obtained by 
division through the length of the cylinders (particles). 
Since the length is rather arbitrary, we drop the unit 
of pressure in the following for the sake of simplicity. 
Due to the linearity of the model, some stress a can be 
rescaled/non-dimensionalized with the stiffness fc2 and 
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particle numbers A'^ 


100, 300 


mean particle radius ao 


1 mm 


relative width of size distribution wq 


±0.5 


particle height h 


6 mm 


particle density pp 


2000 kgm"^' 


particle stiffness k2 


2.105Nm-^ 


particle stiffness k\/k2 


0.5 


particle stiffness kc/k2 


0.5 


particle stiffness fct/fe 


0.2 


particle normal damping 70 


0.2kgs-^ 


particle tangential damping "ft 


0.05 kgs*^ 


particle-particle friction fj. 


0.5 


particle-wall friction /x™ 





melting temperature Tmeit 


120°C 


temperature variation Tvar 


10 K 


density variation with temperature Sut 


10"^ j^-i 


fluid overlap 5**"''' 


0.2 


relaxation time to (default) 


10^ s 


relaxation time during sintering 


10"^ s 


wall mass mw 


0.01 


wall damping 7w 


0.2 



TABLE I: Summary of the system properties and material 
parameters as used in the simulations 



has thus to be read as dimensionless quantity 2.10^a/k2- 
However, a detailed study of the scaling behavior of the 
stresses is far from the scope of this study. 

Note that the pressure is therefore related to the typ- 
ical overlap of two particles: A pressure of p = 100 
corresponds to 100 = 2.10^^2 = 2.10V"/(2aA:2) = 
2.105(5 - So) /{2a) or {S - So)/a « lO^^, while p = 10 
corresponds to approximately {S — So)/a « 10~*. This 
can also be interpreted as mean normal force /" ~ 2ap. 
However, these are rough estimates only, since the con- 
tact forces and the overlaps are strongly varying in mag- 
nitude for one situation, as will be discussed in more 
detail below. 



/ \ ^^'^ 




/ 

/ h 







1 1 
^heat '^sinter 


1 1 
'^cool '^relax 



FIG. 7: Schematic plot of the temperature variation during 
simulation. 

The other system parameters are summarized in table 
H where multiple numbers mean that a series with the 
corresponding values was performed. 



A. Sample preparation 

The preparation of the sample consists of an initial re- 
laxation period at constant temperature T = 80*^0 until 
time ihcat, when the system is heated up to T = 140°C 
between time theat and ^sinter- During the sintering time 
is, the system is allowed to sinter with a much shorter 
relaxation time to = 0.1. This 'trick' allows for a long 
time sintering simulation while keeping the simulation 
time small. During sintering, the time axis should be 
stretched by a factor of 10'' in order to obtain the real- 
time behavior. At the end of the sintering process, at 
time icooi, the sample is slowly cooled down and, at time 
^rciax allowed to relax at constant temperature until time 
tf. With this finished sample, tests will be performed 
later, but first we discuss the preparation process. 

For the preparation of the sample, we use the times 
to = 0, thcat = 0.2, tsintcr - thcat = 0.1, different sintering 
times tg :— tcooi — tsintor, trciax — tcool = 0.1, and tf — 
trciax = 0.1. Note that the time unit is arbitrary, since 
a slower preparation procedure did not lead to noteably 
different results, i.e. the process takes place in the quasi- 
static limit. Only the long-time sintering is affected by a 
change of the relaxation time to- 



1. Temperature and stiffness 

The variation of ki {T, t) with time is plotted in Fig. ||, 
where the algebraic, slow decay of the stiffness with time 
becomes evident. The number of particles was = 300, 
the sidestress p = 100, and the other parameters are 
given in table |. 



2. Density 

The longer the sintering time t^, the lower the value 
of ki gets. For short sintering time, the lowest values 
are never reached, because the system is cooled down 
before the sintering is finished. At the end of the sintering 
time, ki is increasing during the cooling process of the 
sample and reaches its inital value. However, the melting 
and sintering of the contacts is not reversed, as becomes 
evident when plotting the volume fraction i/ in Fig. ^ 

The initial preparation step leads to a rather low den- 
sity of sa 0.8. At time theat, the relaxed sample is heated 
up and, at the same time, the wall friction is switched off. 
The latter has an immediate effect, a slight increase in 
density is possible due to reorganizations. The increased 
temperature becomes only effective at time tsintor, when 
the relaxation time to is decreased in order to acceler- 
ate the evolution of the system. Note that the increase 
in density due to the sintering seemingly appears quite 
rapid - due to the quenched time axis during sintering. 
Up to the beginning of the sintering process all densities 
are equal due to an identical preparation of the sample. 
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FIG. 8: (Top) Temperature during the preparation procedure 
and (bottom) material stifTness ki as function of time for 'ex- 
periments' with different sintering times ts- 



For increasing duration of sintering, ts, the density suc- 
cessivly increases with ts- This effect is strongest at the 
beginning, but continues also for longer times. Below, 
we will discuss different reasons for this increase in den- 
sity in more detail. At the end of the sintering process, 
the system is cooled down, corresponding to the final, 
small increase in density, and then relaxes to the final 
configuration. 



3. Special cases 

As a control simulation, also three particles were simu- 
lation and sintered in the same way as the larger samples 
with N = 100 and N = 300 above. A sample with only 
three particles has almost no possibility for a rearrange- 
ment after the initial load is applied (actually, the three 
particles were either arranged on a triangle or on a line 
and did not change their configuration for the cases ob- 
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FIG. 9: Material density v as function of time for simulations 
with — 300 particles and side stress p — 100 (top) and 
p = 10 (bottom). Note the different vertical axis scaling. 



served). This artificial test-case shows that the strong 
increase in density observed for larger samples is caused 
by re-arrangements of the packing and not only due to the 
force model. More quantitatively, the sintering process 
leads to a densification of about 5% and 16% for p = 10 
and p = 100, respectively, whereas the densification due 
to the contact law, without rearrangements of the pack- 
ing, can be estimated to the respective magnitudes of 
about 0.2% and 8% (data not shown here). Thus, for 
low confining pressure, reorganizations are much more 
important than in the case of large external pressure 
- at least for the parameters used here. In this spirit, 
network-models or DEM simulations with fixed arrange- 
ments ^1, 12, 1^ can not account for the reorganization 
of the packing, which is of eminent importance during 
the sintering process. 

Another control simulation without thermal expansion 
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of the particles, 5aT = 0, gave no new insights and did 
not change the outcome of the simulations after cool- 
down, even though the values of the density were slightly 
different from the values obtained with a thermal expan- 
sion. Therefore, we conclude that the qualitative results 
are not dependent on dctailled parameter values and the 
sintering behavior, macroscopic as well as microscopic, is 
generic within the framework of this model. 



4- Summary 

In summary, we obtain that longer sintering and larger 
confining pressures lead to higher densities of the sintered 
sample. As an example, after the end of the longest sin- 
tering process (dashed line) at i = 1.3, a density of al- 
most 0.92 is reached for p = 100. Only after cool-down, 
the maximum density of almost 0.94 is reached due to 
negative thermal expansion. 

Through test simulations with different parameters, we 
verified that the increase in density is only partially due 
to the contact model, but also is caused by reorganiza- 
tions in the sample. The thermal expansion of the par- 
ticles, on the other hand, does not affect the results in a 
drastic way. 



B. Microscopic picture 

In order to understand the behavior of the material 
during sintering, we take a look at some microscopic 
quantities of the system, like the coordination number, 
in Fig. |l^. 

We note two major issues: First, the larger confining 
stress, p — 100, leads to a larger number of contacts per 
particle, whereas the samples with small confining stress, 
p = 10, approach a coordination number of Cq ~ 4, as 
can be expected for a frictionless, isostatic arrangement 
of disks. Second, the coordination number sometimes 
decays, whereas the density continuously increases. We 
attribute this to the fact, that inside the sample, reor- 
ganizations can take place, which save space and thus 
increase the density, but at the same time can lead to a 
smaller C. 

In order to deepen this insight, the coordination num- 
ber, the density and the measured wall-stress a^z/p are 
plotted in Fig. The continuous increase in density 
is sometimes accompanied by a decrease of the coordi- 
nation number and large stress fluctuations. Large scale 
re-organizations of the packing thus may be responsible 
for a detectable fluctuation of the stress at the boundaries 
of the sample, possibly connected to sound emission. 
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FIG. 10: Coordination number as function of time during the 
preparation and sintering {t,, = 1 s) of the sample. Compared 
are samples of different size, = 100 (a) and TV = 300 (b), 
and with different confining stress p. 



C. Contact statistics 

1. Compression/deformation/overlap probabilities 

A very specific microscopic property of the sintered 
sample is the statistics of the contacts, i.e. the probability 
distribution function (pdf) to find a certain overlap S. 
After the initial preparation of the powder sample at low 
temperature, the probability for larger overlaps decays 
rapidly - no large deformations have taken place. 

The overlap probabilities at the final stage of the sin- 
tering process are plotted in Fig. 
both side pressures p = 10 and p = 100. 



12| for TV = 300, and 
The probability 
distributions show many contacts with small overlaps and 
a rapidly decaying probability for large overlaps at short 
sintering times. The mean overlap increases, as expected. 
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FIG. 11: Density (squares), coordination number (circles), 
and normalized stress on the top- wall (line), as function of 
time during the sintering {ts = 1 s) of the sample. The coor- 
dination number is scaled by a factor 0.2 in order to allow a 
comparison with the other quantities. 



with the side stress and the sintering time. The larger 
side stress leads to an underpopulation of small overlaps, 
because all particles are compressed quite strongly. 

The probability distributions after the cool-down and 
the subsequent relaxation is plotted in Fig. |l^ for the 
same simulations. The remarkable difference is not the 
width of the distribution - that is only slightly wider. 
The most striking difference is the shape, showing that 
after cool-down contacts with small overlap are rarefied, 
while contacts with a the typical overlap are overpop- 
ulated as compared to the situation before cool-down. 
The probability for the largest overlaps is still decreasing 
rapidly. 



2. Normal contact force probabilities 

Due to the advanced contact force law, as introduced 
and used in this stude, it is not possible to directly obtain 
the contact force from the overlap. Therefore, the pdf for 
the normal contact forces is plotted in Figs. |lj and |l^ 
for small and large side-stresses p = 10 and p = 100, 
respectively, and for different sintering times, tg as given 
in the inset. 

Comparing the pdf for the two situations, after sinter- 
ing and after cool-down, the astonishing outcome of the 
simulation is the fact that the contact forces are mostly 
repulsive and rather small in the hot situation, just be- 
fore cool-down. This situation changes during cooling 
down: The cooling goes ahead with a broadening of the 
distribution towards both positive and negative forces. 
This qualitative result is independent of the side pressure, 
however, the width of the distribution increases with in- 
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FIG. 12: Overlap probability for overlaps, rescaled with the 
sum of the radii ai and a2 of the contacting particles, after 
different sintering durations. The data are taken from the 
simulations with N = 300 at different confining pressure p. 



creasing side stress and increasing sintering time. 

Another interesing and unexpected result is the obser- 
vation that the force distribution becomes narrower with 
sintering. In other words, the extremely large forces are 
"destroyed" due to long time sintering and the distri- 
bution becomes more homogeneous in the sense that its 
width becomes smaller, see top panel in Fig. 15 



3. Tangential forces 

From the tangential force distribution (data not 
shown), there are less clear observations to be made. 
The force distribution shows that during sintering the 
tangential forces become weaker - as can be expected 
from the force model, see Eq. (24), since ki{T,t) decays 
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FIG. 14: Contact force probability distribution after the sin- 
tering (top) and after the cool-down and relaxation (bottom), 
for iV = 300 and p = 10. 



FIG. 13: Overlap probability for overlaps, rescaled with the 
sum of the radii ai and 02 of the contacting particles, after 
different sintering durations. The data are taken from the 
simulations with A'^ = 300 at different confining pressure p. 



with increasing temperature and time. Thus longer sin- 
tering leads to weaker tangential forces. However, after 
cool-down, the tangential force distribution is almost in- 
dependent of the side stress and the sintering time, if 
both are sufficiently large. With other words, only for 
very small side stresses and sintering times, the tangen- 
tial forces are weakly activated, for larger p or t^, the 
tangential forces reach a saturation distribution that is 
decaying almost exponentially for large ft- 

The second possibility to look at the tangential forces 
is to measure the amount of friction that is activated, 
namely fit := ft/ifn - /min)- This quantity should not 
become larger than ii{T,t), as is consistently observed 
from the data. Contacts where fit ~ fJ-iT-it) are referred 
to as contacts with fully activated friction. These be- 
come less probable for longer sintering times since, as 



mentioned before, the attractive forces become stronger 
after long sintering and the cooling down. Larger attrac- 
tive forces correspond to a larger magnitude of /min, so 
that fit becomes smaller. 



D. Material properties 

The samples prepared via the procedure described 
above are now tested via two test experiments, a com- 
pression test, where the vertical confining pressure is 
slowly increased, and a vibration test, where the con- 
fining stress is removed, and the sample is vibrated ver- 
tically on a flat bottom. The compression test is per- 
formed for two sample sizes, namely with = 100 and 
iV = 300, in order to judge the effect of a rather small 
particle number on the outcome of the simulations. 
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higher compressive pressure leads to higher densities and 
also to a broader variation in densities between. 
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FIG. 15: Contact force probability distribution after the sin- 
tering (top) and after the cool-down and relaxation (bottom), 
for iV = 300 and p = 100. 




1. Compression test (100 particles) 

The compression test is a variant of the bi-axial com- 
pression frequently used in soil- or powder mechanics. A 
given compressive stress in one direction {p — axx) is kept 
constant by moving the right, vertical wall, if neccessary. 
The other confining stress a^z is increased by moving the 
vertical wall down in a defined way - see also section ||. 
This vertical motion is thus strain controlled, where the 
vertical strain is defined as tzz{t) = 1 — z{t)/ Zf), with the 
vertical position of the top wall z{t) and zq = z{Q). Note 
that a compressive vertical strain is defined positive for 
convenience here. 

For the compression test, the previously prepared sam- 
ples are used and the top wall is displaced slowly with 
the vertical strain ezz{t)- The density of a sample with 
A'' = 100 particles is plotted against the strain for vari- 
ous sintering times and for the two pressures p = 10 and 
p = 100 in Fig. |l6|. The results show again that samples, 
which sintered longer have a higher density. Further- 
more, the compressive pressure also affects the sample 
density in magnitude and also in the variation, i.e. the 
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FIG. 16: Material density as function of the vertical strain 
for samples prepared with different sintering times ts and dif- 
ferent pressure pw = 10 (Top) and pw = 100 (Bottom). 

During compression, the density slightly increases first 
and then decreases strongly. The former is due to com- 
pression, the latter is due to dilatancy that is neccessary 
for the material to shear. The longer the sample was sin- 
tered, the stronger is the change in density. We relate 
this to the existence of the attractive contact forces. For 
longer sintering, stronger attraction is activated, holding 
together parts of the sample, so that the compression test 
leads to fragements with increasing size for increasing sin- 
tering duration. For the shortest sintering, the density 
change due to compression is negligible and the sample 
fragments into single grains. 

The behavior of the material, as for instance its stiff- 
ness, is better described in terms of the vertical stress 
that the material can sustain under load, as plotted in 



14 



Fig. The vertical stress in the sample increases and 
it fails typically at some strain and at some magnitude 
of stress. The failure stress increases with increasing sin- 
tering time and increasing external stress. The material 
stiffness (the slope in this representation) is increased by 
a factor of about two when the confining stress is in- 
creased by a factor of ten. Moreover, the critical strain 
where the material fails increases with increasing con- 
fining pressure. Finally, a rather large jump in material 
strength is observed for sintering times between ts = 0.02 
and 0.05 for the small pressure, not paralleled by a simi- 
lar outcome for the large stress. Thus the combination of 
sintering time, compressive pressure and test mode ap- 
pears quite non-linear and not straightforward. 
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2. Compression test (300 particles) 

For the compression test of the larger samples, also the 
prepared samples are used and the top wall is displaced 
slowly with the vertical strain. The density of a sample 
with N = 300 particles is plotted against the strain for 
various sintering times and for the two pressures p in 
Fig. |l^. The results are qualitatively similar, only that 
the larger sample shows somewhat larger densities and a 
weaker decay of density after failure. 
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FIG. 18: Material density as function of the vertical strain 
for samples prepared with different sintering times ta and dif- 
ferent pressure pw = 10 (Top) and pw = 100 (Bottom). 



FIG. 17: Vertical stress as function of the vertical strain for 
samples prepared with different sintering times ts and differ- 
ent pressure pw = 10 (Top) and pw = 100 (Bottom). The 
slope of the dashed line (material stiffness) is about double 
for the higher pressure. 



The vertical stress is plotted against the strain in Fig. 
[l9| . With increasing strain, the vertical stress in the sam- 
ple increases and it fails typically at some slightly larger 
strain and stress as in the small sample. The failure 
stress increases with increasing sintering time, increasing 
external stress, and increasing sample size. However, the 
latter observation is not clearly a size-effect, since the 
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sample size is so small that a non-negligible fraction of 
the particles is in contact with the walls and thus leads 
to different outcome. The material stiffness (dashed line 
in Fig. |l^) is again increased by about a factor of two 
when the confining stress is increased by a factor of ten. 
Moreover, the critical strain where the material fails in- 
creases with increasing confining pressure, sintering time 
and system size. 
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FIG. 19: Vertical stress as function of the vertical strain for 
samples prepared with different sintering times ts and differ- 
ent pressure pw = 10 (Top) and Pw = 100 (Bottom). The 
slope of the dashed line (material stiffness) is about double 
for the higher pressure. 



3. Compression test snapshots 

One compression test for a long sintering time, N = 
300 and p = 10 is presented in Fig. ^ During compres- 
sion (from top to bottom), the lines of strong attraction 



(so to say the backbone of the sample) are desctroyed 
and. at the same time, more and more frictional con- 
tacts occur due to local shear. Moreover, gaps between 
the parts of the sample open and it fragments into pieces. 

This is only a representative example for the compres- 
sion test; a more detailed study of the fracture behavior, 
sample-size- and sintering-time-dependence is far from 
the scope of this study. 




FIG. 20: Snapshots of the compression test with = 300 and 
p = 10. at times t = 0.005, 0.02, 0.03, 0.035, and 0.04 (from 
top to bottom); the material sample was sintered for ts = 1.0. 
The blue bars are the walls, the circles are the particles with 
their colors coding the average stress, blue, green and red cor- 
respond to low, medium and large stresses. The lines are the 
contacts with their colors coding attractive (blue) or repulsive 
(red) normal forces. The small solid circles denote the tan- 
gential forces, with their size proportional to the magnitude 
of the tangential force. 
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E. Vibration test 

The sintered samples can also be vibrated (in the grav- 
itational field) in order to probe their stability. The re- 
sults of the previous tests are paralleled by a vibration 
test, see Figs. ^ and |2^. 



pies, whereas the sample becomes more and more sta- 
ble with increasing sintering time and confining pressure. 
The sample with the longest sintering time t<; = 1.0 is al- 
most perfectly stable even under strong shaking - some 
corner- or boundary particles sometimes break off. For 
shorter sintering time the sample is less stable and frag- 
ments into pieces. For very short sintering time, the sam- 
ple consists of single particles only. 
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FIG. 21: Snapshots after a vibration test of t — 0.1, with fre- 
quency / — 100 Hz and amplitude a = 0.2 mm. The material 
sample was sintered (from top to bottom) for ts — 10"^, 0.02, 
and 1.0 with p — 10. The open circles are particles with their 
colors coding the average stress, blue, green and red corre- 
spond to low, medium and large stresses. The lines are the 
contacts with their colors coding attractive (blue) or repulsive 
(red) normal forces. The small solid circles denote the tan- 
gential forces, with their size proportional to the magnitude 
of the tangential force. 



FIG. 22: Snapshots after a vibration test of t = 0.1, with fre- 
quency / = 100 Hz and amplitude a = 0.2 mm. The material 
sample was sintered (from top to bottom) for ts — 10~*, 0.01, 
and 1.0 with p = 100. The color coding and figure meaning 
is the same as in Fig. 



Short sintering times leads to unstable material sam- 
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V. CONCLUSION 

In summary, a discrete model for the sintering of par- 
ticiilatc materials was introduced and simple material 
samples were sintered for different times and confining 
pressures. Then they were tested with respect to their 
anisotropic load strength: Longer sintering and stronger 
confining pressure systematically increases the density 
and the strength of the material. Depending on the sin- 
tering duration, either isolated particles, fragments or a 
single solid block of material could be produced. 

Besides this macroscopic point of view, also the mi- 
croscopic picture was examined. The series of astonish- 
ing observations includes: (i) The coordination number 
may slightly decrease due to reorganizations while the 
density monotonously increases, (ii) Sintering leads to 
a broadening overlap distribution, but to a narrowing 
force distributon, and thus to a homogenization of the 
sample, (iii) The normal forces become strongly attrac- 
tive during the cooling down of the sample below the 
melting temperature. Besides these facts, a lot of open 
questions concerning the sintering process remain, espe- 
cially concerning the connection between the microscopic 
contact-model parameters and the macroscopic material 
parameters. 



The research to be done is an accurate testing of the 
model via a comparison with experimental data. Since 
these are only available in three-dimensional systems, the 
2D model presented here may be not helpful. However, 
the model can easily be extended to three dimensions, 
where 'only' more particles are needed. Note that the 
model presented here increases the amount of compu- 
tation necessary for each contact by a large factor, so 
that the number of particles possible to simulate becomes 
rather small for a standard computer. Thus an extension 
to three dimensions requires a proper tuning of the im- 
plementation of the force-model and possibly the view 
of particles as blocks of material that cannot fragment, 
rather than isolated particles. The last missing ingredi- 
ent in the model is a rolling resistance which accounts for 
a torque resistance of the contacts. 
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